Correction of the effect of gradient field non-linearities in phase contrast MRI

ABSTRACT

Errors in qualitative phase contrast measurements due to gradient field heterogeneities are reduced by using either a generalized reconstruction algorithm or an approximate reconstruction algorithm. True velocities are calculated using measured velocity information and phase differences, first moments of gradients, and gyromagnetic ratio.

CROSS-REFERENCE TO RELATED APPLICATION

[0001] This patent application is a continuation in part of co-pendingU.S. patent application No. 10/317,516 filed Dec. 11, 2002, for“Correction of the Effects of Spatial Gradient Field Distortions inDiffusion-Weighted Imaging,” (Atty. Dkt. No. STFUP096/S02-166), which isincorporated herein by reference for all purposes.

GOVERNMENT RIGHTS

[0002] The U.S. Government has rights in this invention pursuant to NIHGrant No. P41RRO9784 to Stanford University.

BACKGROUND

[0003] This invention relates generally to magnetic resonance imaging(MRI), and more particularly the invention relates to flow imaging usingphase contrast MRI.

[0004] Magnetic resonance imaging (MRI) requires placing an object to beimaged in a static magnetic field, exciting nuclear spins in the objectwithin the magnetic field, and then detecting signals emitted by theexcited spins as they precess within the magnetic field. Through the useof magnetic gradient and phase encoding of the excited magnetization,detected signals can be spatially localized in three dimensions.

[0005] Phase contrast MRI (PC-MRI) is widely used to assess blood flowand tissue motion. The technique relies on the measurement of changes inthe signal phase due to flow or motion in the presence of known linearmagnetic gradient fields.

[0006] Although it is well known that non-uniformity in magnetic fieldgradients can cause significant image warping and require correction,little has been reported about the impact of spatial gradient fielddistortions on velocity encoding.

[0007] In PC-MRI, these imperfections introduce errors in velocitymeasurements by affecting the first moments used to encode flow ormotion. Typically, based on the deviations of the actual gradient fromthe desired, uniform gradient, the spatial image distortions inmagnitude and phase images are retrospectively corrected by imageunwarping algorithms, e.g., with algorithms included in the imagereconstruction software. The velocity-encoded information is moved toits correct location but the error in velocity encoding due to the localfield deviation still persists. Previous work (e.g., U.S. Pat. No.6,163,152) recognized the impact of errors in the strength of theencoding gradient and disclosed a correction for this effect. In somerecent high performance gradient systems, the coil size was reduced tolimit dB/dt and amplifier power. As a result, gradient uniformity hasbecome even more degraded.

SUMMARY OF THE INVENTION

[0008] The inventors recognized that error in not only the strength butalso the direction of the local gradient from its ideal value isdirectly reflected in errors in the first order gradient moments andthus the velocity encoding. Therefore, gradient field distortions canlead to c not only deviations from the nominal gradient strength butalso from the intended gradient direction and thus affects not only themagnitude of encoded velocities but also velocity encoding direction inPC MRI.

[0009] In accordance with the invention, a generalized method for thecharacterization of gradient field non-uniformity and accuratereconstructions of phase contrast data is presented. While acquisitionof three-directional velocity information is required to apply thisgeneralized reconstruction for typical applications of phase contrastMRI, the full three-directional velocity information may not be readilyavailable or difficult to obtain for practical reasons (e.g., limitedtotal acquisition time in breath-held measurements or desired hightemporal resolution for a single velocity component such as throughplane flow at the level of heart valves). Therefore, an additionalapproximate solution which does not require the full three-directionalvelocity information is provided for practical applications.

[0010] The invention and objects and features thereof will be morereadily apparent from the following detailed description and appendedclaims when taken with the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

[0011] FIGS. 1A-1C are simulated iso-contours (labeled in units of mT)for the magnetic field B_(z,model) ^((i)) (x, z) generated by thei-gradient coil (i=x,y,z, FIG. 1A: z-gradient coil, FIG. 1B: x-gradientcoil, FIG. 1C: y-gradient coil) in coronal planes at y=0 mm and y=200 mm(left column) with the overlaid vector field representing the projectionof the magnetic field gradient into the coronal plane.

[0012]FIGS. 2A, 2B show results of 3D phase contrast MRI measurementswith three-directional velocity encoding including two dimensionalintensity plots reflecting velocity profiles in the coronal (x-z) planefor FIG. 2A) velocities encoded along z, and FIG. 2B) velocities encodedalong y.

[0013] FIGS. 3A-3C show results of 3D phase contrast MRI measurementswith three-directional velocity encoding including measured velocitiesv_(x), v_(y) and v_(z) along cross sections through the center of thetube phantom (as indicated by the dashed line in the right plot in FIG.2B) with each graph showing velocities for standard (solid gray line,standard PC) and generalized phase contrast reconstruction (thick blackdashed line, generalized PC) as well as the difference between bothmethods (black dot-dashed line).

[0014]FIGS. 4A, 4B show results of 2D phase contrast MRI measurementswith one-directional velocity encoding and demonstrates the effect ofgradient field inhomogeneities on velocity encoding for locations of thetube phantom as defined by the (x,y) coordinates.

[0015]FIG. 5 illustrates 2D phase contrast MRI with one-directionalvelocity encoding including correlation of measured versus model baseddeviations Δ{overscore (v)}_(z) from through plane velocities.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0016] As noted above, inhomogeneous, or heterogeneous, gradient fieldscan introduce deviations from the nominal gradient strength andorientation and therefore spatially dependent first gradient moments.Resulting errors in the measured phase shifts used for velocity encodingcan therefore cause significant deviations in velocity quantification.In accordance with the invention, the true magnitude and direction ofthe underlying velocities can be recovered from the phase differenceimages by a generalized phase contrast velocity reconstruction whichrequires the measurement of full three-directional velocity information.The generalized reconstruction of velocities is applied using a matrixformalism that includes relative gradient field deviations derived froma mathematical model of local gradient field non-uniformity. Inaddition, an approximate solution for the correction of one-directionalvelocity encoding is provided. Dependent upon the spatial location ofthe velocity measurements, the magnitude of encoding errors can be ashigh as 60% while errors in the encoding direction can be up to 45°.Results of phantom measurements demonstrate that effects of gradientfield non-uniformity on PC-MRI can be corrected in accordance with theinvention.

[0017] To model the spatial imperfections of the magnetic fieldB_(z,model) ^((i))(r) produced by the gradient coils along nominaldirections i=x, y, and z, a polynomial model based on a sphericalharmonic expansion is generally employed [6], although other models canalso be used. Scanner specific parameters are included in thecoefficients used to describe the spatial dependence of the gradientfield (appendix, equation A1).

[0018] Based on the model, local magnetic field gradientsG_(ji,model)(r) can be calculated according to $\begin{matrix}{{{G_{{ji},{model}}(r)} = {\frac{\partial{B_{z,{model}}^{(i)}(r)}}{\partial j}\quad i}},{j = x},y,{{and}\quad {z.}}} & (1)\end{matrix}$

[0019] In this context G_(ji,model)(r) describes the component of thegradient along j for an ideal (nominal) gradient along i.G_(ji,model)(r) therefore is the component of the actual gradient fieldalong the selected gradient field direction (self term), whileG_(ji,model)(r) with i≠j are gradient terms orthogonal to the primarilyselected gradient direction (cross terms).

[0020] As a result, the true gradient field demonstrates not onlydeviations from the nominal gradient strength but also from the originalgradient direction i. The relative deviation from the ideal (spatiallyconstant) gradient field G_(i,ideal) of the actual (modeled) fieldstrength G_(ji,model)(r) is: $\begin{matrix}{{{\Lambda_{ji}(r)} = {\frac{G_{{ji},{model}}(r)}{G_{i,{ideal}}}\quad i}},{j = x},y,{{and}\quad z}} & (2)\end{matrix}$

[0021] for each physical gradient axis i=(x, y, z). Here, the vectorr=(x,y,z)^(T) reflects the true (already spatially unwarped) position.

Generalized Reconstruction

[0022] The relative field deviations Λ_(ji)(r) result in a relativeerror in the associated first gradient moments and therefore velocityinduced phase shift in phase contrast MRI. With the ideal first momentM_(1,ideal) for an ideal arbitrary gradient G_(ideal)(t) given by$\begin{matrix}{M_{1,{ideal}} = {\int_{0}^{TE}{{G_{ideal}(\tau)}\tau {\tau}}}} & (3)\end{matrix}$

[0023] and equations (1) and (2) the actual first moment M₁(r) is nowspatially dependent and can be calculated using the matrix Λ(r) whichcontains the relative field deviations Λ_(ji)(r): $\begin{matrix}\begin{matrix}{{M_{1}(r)} = {\begin{pmatrix}{M_{1,x}(r)} \\{M_{1,y}(r)} \\{M_{1,z}(r)}\end{pmatrix} = {\int_{0}^{TE}{\begin{pmatrix}{\Lambda_{xx}(r)} & {\Lambda_{xy}(r)} & {\Lambda_{xz}(r)} \\{\Lambda_{yx}(r)} & {\Lambda_{yy}(r)} & {\Lambda_{yz}(r)} \\{\Lambda_{zx}(r)} & {\Lambda_{zy}(r)} & {\Lambda_{zz}(r)}\end{pmatrix}\begin{pmatrix}{G_{x,{ideal}}(\tau)} \\{G_{y,{ideal}}(\tau)} \\{G_{z,{ideal}}(\tau)}\end{pmatrix}\tau {\tau}}}}} \\{= {{{\Lambda (r)}{\int_{0}^{TE}{{G_{ideal}(\tau)}\tau {\tau}}}} = {{\Lambda (r)}M_{1,{ideal}}}}}\end{matrix} & (4)\end{matrix}$

[0024] For a typical phase contrast experiment, the actually measuredphase difference Δφ(r) with encoding along a single (arbitrary)direction (assuming constant velocity v(r) and a gradient waveform withzero net area at echo time TE) is therefore given by $\begin{matrix}\begin{matrix}{{{\Delta\varphi}(r)} = {{v(r)}^{T}{\gamma \left\lbrack {{{\Lambda (r)}M_{1,{ideal}}^{(I)}} - {{\Lambda (r)}M_{1,{ideal}}^{({II})}}} \right\rbrack}}} \\{= {{v(r)}^{T}{{\gamma\Lambda}(r)}\Delta \quad M_{1,{ideal}}}} \\{= {{\gamma \left\lbrack {{\Lambda (r)}\Delta \quad M_{1,{ideal}}} \right\rbrack}^{T}{{v(r)}.}}}\end{matrix} & (5)\end{matrix}$

[0025] The two different ideal first moments along one gradientdirection which are necessary to encode a single flow direction aredenoted by M_(1,ideal) ^((I)) and M_(1,ideal) ^((II)); γ is thegyromagnetic ratio. The velocity sensitivity of the experiment isdetermined by the difference ΔM_(1,ideal) in first moments.

[0026] For the most general case, a phase difference vectorΔΦ(r)=(Δφ₁(r),Δφ₂(r), . . . , Δφ_(N)(r))^(T) resulting from velocityencoding along N different (arbitrary) directions (with nominal velocitysensitivities characterized by an (N×3) encoding matrix Ω=(ΔM_(1,ideal)⁽¹⁾,ΔM_(1,ideal) ⁽²⁾, . . . , ΔM_(1,ideal) ^((N)))) can be written as

ΔΦ(r)=γ[Λ(r)Ω]^(T) v(r).  (6)

[0027] Note that in this notation, superscripts on M₁ (e.g., M_(1,ideal)^((II)), M_(1,ideal) ^((II)) identify first moments of subsequentindividual measurements, whereas superscripts on ΔM₁ (e.g., ΔM_(1,ideal)⁽¹⁾,ΔM_(1,ideal) ^((N))) denote velocity sensitivities (differences infirst moments) corresponding to different encoding directions.

[0028] For a general correction for gradient field inhomogeneities thetrue velocities of underlying flow or motion can therefore be recoveredfrom the measured phase differences by inverting [Λ(r)Ω]^(T) and solvingequation (6) for v(r): $\begin{matrix}{{v(r)} = {{\frac{1}{\gamma}\left\lbrack {{\Lambda (r)}\Omega} \right\rbrack}^{- T}\Delta \quad {{\Phi (r)}.}}} & (7)\end{matrix}$

[0029] In the general case, the true velocity vector can only berecovered from the phase difference data if at least three-directionalvelocity encoding (N=3) was performed. A more general discussion of thesolutions of equation (6) and the role of Singular Value decomposition(SVD) in the over- or underdetermined case is given in the appendix.

[0030] In conventional reconstruction of velocity data from phasecontrast measurements, on the other hand, gradient field deviations areignored and velocity components are calculated simply by scaling thephase difference images according to the magnitude of the ideal firstmoment differences ΔM_(1,ideal) ^((i)) for each encoding directionindependently: $\begin{matrix}{{{v_{i}(r)} = {\frac{1}{{\gamma\Delta}\quad M_{1,{ideal}}^{(i)}}\Delta \quad {\varphi_{i}(r)}}},{i = 1},\ldots \quad,{N.}} & (8)\end{matrix}$

[0031] A more detailed description of the standard reconstruction isgiven in Bernstein et al., “Reconstructions of Phase Contrast, PhasedArray Multicoil Data,” Magn. Reson. Med. 1994; 32(3):330-4.

[0032] Implications of the generalized reconstruction (equation (7)) ontypical applications of phase contrast MRI such as velocity encodingalong a single axis (N=1) and along three orthogonal axes (N=3) aregiven in the appendix.

Single-Direction Encoding

[0033] In some applications (e.g., single-direction encoding which isquite frequently used to measure volumetric flow rate) the fullthree-directional velocity information is not available (see alsoequation A2, appendix). In these cases, one is usually interested in thethrough plane component of the velocity in order to determine the flowalong this direction. Due to the gradient non-uniformities, the flowencoding is changed in both magnitude and direction. However, the localslice select direction is also altered and is parallel to the actual,local, flow encoding direction. This parallel relationship is ideal forflow quantitation, and for this application, it is sufficient to correctfor the magnitude error in the flow encoding. Suppose the intended flowencoding (and slice select) direction is along the direction of the unitvector n with the ideal first moment for single-direction encodingrepresented by ΔM_(1,ideal) ⁽¹⁾=ΔM_(1,ideal) ⁽¹⁾n. In general, theactual velocity encoding and slice select direction is given byn′(r)=Λ(r)n/∥Λ(r)n∥ and the velocity component normal to the excitedslice at the location r along the actual velocity encoding direction isv_(⊥)(r)=n′(r)^(T)v(r)=[Λ(r)n]^(T)v(r)/∥Λ(r)n∥. For single-directionvelocity encoding equation (6) becomes

Δφ₁(r)=γ[Λ(r)ΔM _(1,ideal) ⁽¹⁾]^(T) v(r)=γΔM _(1,ideal) ⁽¹⁾[Λ(r)n] ^(T)v(r)=γΔM _(1,ideal) ⁽¹⁾ v _(⊥)(r)∥Λ(r)n∥  (9)

[0034] and the velocity normal to the local excited slice is$\begin{matrix}{{v_{\bot}(r)} = \frac{{\Delta\varphi}_{1}(r)}{{\gamma\Delta}\quad M_{1,{ideal}}^{(1)}\sqrt{\sum\limits_{j = 1}^{3}\left\lbrack {\sum\limits_{i = 1}^{3}{n_{i}{\Lambda_{ji}(r)}}} \right\rbrack^{2}}}} & (10)\end{matrix}$

[0035] For the special case where n points along i=x,y, or z, equation(10) reduces to $\begin{matrix}{{v_{\bot}(r)} = \frac{{\Delta\varphi}_{1}(r)}{\gamma \quad d_{i}\Delta \quad M_{1,{ideal}}^{(1)}}} & (11)\end{matrix}$

[0036] with d_(i)=(Λ_(ii) ²+Λ_(ji) ²+Λ_(ki) ²)^(1/2) where j and k arethe corresponding orthogonal directions. Note that equations (10) and(11) offer accurate reconstruction of v_(⊥)(r), but the actual flowdirection and speed are generally not known.

[0037] In other applications of single-direction encoding the velocitycomponent along the intended encoding direction n (v_(n)(r)) may be ofinterest. It is important to note that, in general, it is impossible toexactly reconstruct v_(n)(r) in the presence of gradient fieldinhomogeneity. However, an approximate solution may be useful in somecases when the direction error is small or if it can be assumed that thetrue flow is predominantly along n (i.e., v(r)=v_(n)(r)n). In thisapproximation the velocities v_(n)(r) can be estimated by$\begin{matrix}{{v_{n}(r)} = \frac{{\Delta\varphi}_{1}(r)}{{\gamma\Delta}\quad M_{1,{ideal}}^{(1)}{\sum\limits_{i,j}^{\quad}{n_{i}n_{j}\Lambda_{ij}}}}} & (12)\end{matrix}$

[0038] Note that, for encoding along one of the Cartesian directionsi=x,y or z the sum in equation (12) reduces to the self term Λ_(ii)(r)along the intended encoding direction i.

[0039] The case described above pertained to the application in which asingle velocity component was desired and no a-priori knowledge of thetrue flow direction was available. In the event the true flow directionis known or some assumptions about velocity ratios can be made (e.g.,v_(z)=mv_(y) and v_(z)=lv_(x)), the cross terms can be included in thereconstruction and one-directional velocity encoding along the intendedencoding direction n can be fully corrected using equation (7).

[0040] In demonstrating the invention, all measurements were performedon a 1.5T system (Signa CV/i, GE, Milwaukee, Wis., G_(max)=40 mT/m,risetime=268 μs).

[0041] The phantom used for all measurements consisted of a tube(diameter=20 mm), parallel to the direction of the main magnetic field(z), which could be moved to different (x-y) locations in the axialplane. A computer controlled pump was used to generate a constantlaminar flow at a rate of 25 ml/s. A long tube (2 meters) was used toensure that velocity profiles were fully developed before they enteredthe imaging volume.

[0042] In order to demonstrate velocity encoding errors related togradient field non-uniformity and the feasibility of both thegeneralized reconstruction (equation (7)) and the approximate solution(equation (12)) two sets of phantom experiments were performed with fullthree-directional and one-directional through-plane velocity encoding,respectively.

[0043] For all experiments, the 3D volume or a set of 2D slices wereprescribed within a volume so as to cover a large extent of gradientfield inhomogeneities. Data acquisition was performed with an in-planeimage matrix of 256×256, a FOV of 40-cm and velocity encoding (venc) of30 cm/s.

[0044] For the analysis of the generalized phase contrast reconstruction3D scans were performed with a coronal (x-z) slab using a RF-spoiledgradient echo sequence (TE=3 ms, TR=10.3 ms, BW=+/−62.5 kHz, α=15°) withidentical velocity encoding (venc=30 cm/s) along all three spatialdirections. The slab thickness was adjusted to include the whole tube(slab thickness=72 mm, 24 slice encodes).

[0045] During a second set of phantom measurements a 2D RF-spoiledgradient echo sequence with velocity encoding along a single directionwas used (TE=3.9 ms, TR=9 ms, α=20°) to demonstrate the feasibly of theapproximate correction (equation (12)). Twenty-one equally spaced axialimages (slice thickness=10 mm, gap=10 mm) were acquired withthrough-plane flow encoding only.

[0046] In order to correct for the effect of eddy current and systemimperfections other than gradient inhomogeneities, the measured phasedifference images for each experiment was referenced to a scan with thepump off (i.e., no flow) but otherwise identical parameters.

[0047] All gradient field calculations and data reconstruction wereperformed using MATLAB (The MathWorks, USA) on a personal computer.Images were computed with the proposed generalized and also withstandard phase contrast reconstruction algorithms.

[0048] The magnitude and velocity images from the 3D data set werecompletely spatially unwarped (in 3D). The data were then used togenerate two-dimensional intensity plots and cross sectional velocityprofiles for generalized and standard phase contrast reconstruction.

[0049] For the second set of phantom experiments, images were unwarpedin 2D and the data were quantitatively analyzed by calculating the meanthrough plane velocities for the 21 equally spaced axial slice planes.Velocities were calculated within ROIs that were defined at thedifferent tube locations in the x-y-plane for each slice. Since gradientfield distortions also affect slice plane locations and are notcorrected by the scanner's unwarping software (the correction is onlydone in the in-plane dimensions) the slice (z) position for eachcalculated {overscore (v)}_(z) for each tube was also corrected formis-mapping and moved back to its true location.

Gradient Field Distortions

[0050]FIGS. 1A, 1B, and 1C are simulated iso-contours (labeled in unitsof mT) for the magnetic field B_(z,model) ^((i))(x, z) generated by thei-gradient coil (i=x,y,z, A: z-gradient coil, B: x-gradient coil, C:y-gradient coil) in coronal planes at y=0 mm and y=200 mm (left column).The overlaid vector field represents the projection of the magneticfield gradient into the coronal plane. Angular deviations θ_(i)(x,z) indegrees and relative strength variation d_(i)(x,z) from the nominalgradient field are displayed as gray-scaled images. All plots reflectgradient fields within coronal planes with a FOV of 400 mm; x, y and zcoordinates are given in mm. Note that, according to the model, gradientfield distortions for gradient coil axes along x and y axis areidentical but transposed.

[0051]FIG. 1 illustrates the spatial dependence of gradient fielddistortions derived from a gradient field model given including thecalibration coefficients of our system. If the z-gradient coil isenergized the spatial gradient field imperfection can be characterizedby the local angular deviation from the nominal gradient directionθ_(z)=tan⁻¹((Λ_(xz) ²+Λ_(yz) ²)^(1/2)/Λ_(zz)) and the local absolutedeviation from nominal gradient strength d_(z)=(Λ_(zz) ²+Λ_(xz) ²+Λ_(yz)²)^(1/2). Analogous expressions are defined for x- and y-coils.

[0052] FIGS. 1A-1C show iso-contours of the modeled magnetic fieldB_(z,model) ^((i))(x, z) generated by the i-gradient coil (i=x,y,z) incoronal planes through the iso-center of the magnet (y=0 mm, top row)and off-center at y=100 mm (bottom row). The overlaid vector fieldrepresents the projection of the magnetic field gradient into thecoronal plane. Angular and absolute deviations from the nominal gradientfield are illustrated by intensity plots of θ_(i)(x,z) and d_(i)(x,z),respectively. Note that FIG. 1C depicts gradient fields with a principalgradient direction perpendicular to the coronal plane.

[0053] Ideal gradient fields would result in equidistant contour linesof the magnetic field in FIGS. 1a and 1 b and spatially constantB_(z,model) ^((y))(x, z) in FIG. 1C. Deviations from those idealconditions are clearly visible for all three gradient axes. Inprinciple, based on the gradient field model, encoding errors can be ashigh as +/−60% in magnitude and up to 45° in angular deviations within avolume of (400 mm)³, depending on the spatial location and encodingdirection. Table 1 lists maximum angular and absolute deviations for allthree gradient coils (x and y coils are identical but rotated) withintwo different volumes ((200 mm)³ and (400 mm)³) centered at theiso-center of the magnet. Note that the deviations from nominal gradientstrength d_(i) directly correspond to locally varying errors inthrough-plane velocity components v_(⊥)(r) for single-directionencoding, which can be fully corrected using equation (11). TABLE 1Maximum deviations in gradient strength d_(i) and orientation θ_(i) (i =x, y, z) for measurement volumes of (200 mm)³ and (400 mm)³. max(θ_(i))[deg] max(d_(i)) [%] x, y-gradient z-gradient x, y-gradient z-gradientcoil coil coil coil (200 mm)³ 13.2° 1.1° 12% 4% volume (400 mm)³ 45.1°16.2° 63% 67%  volume

[0054] From both FIG. 1 and Table 1 it is evident that deviations in thestrength of the encoding typically increase with distance from theiso-center of the magnet while the angular deviation demonstrates morecomplex patterns. In addition, the x- and y-gradient coils producestronger deviations from the nominal fields than the z-gradient coil.Symmetries in the gradient field model are a result of the fact that thephysical x- and y-gradient coils are identical but rotated.

Generalized Reconstruction

[0055]FIGS. 2A and 2B illustrate 3D phase contrast MRI withthree-directional velocity encoding: two dimensional intensity plotsreflecting velocity profiles in the coronal (x-z) plane for A)velocities encoded along z and B) velocities encoded along y. The centerof the tube phantom was close to the iso-center of the magnet (y=−20mm). Each set of images shows velocity profiles for standard (left,standard PC, equation (8)) and generalized phase contrast reconstruction(middle, generalized PC, equation (7)) as well as the difference betweenboth methods (right).

[0056]FIG. 2 depicts results of the 3D phantom measurement withthree-directional velocity encoding. The intensity plots show pixel wisevelocity profiles calculated by standard (equation (8)) and generalizedphase contrast reconstruction (equation (7)) as well as differenceimages. All images correspond to a thin coronal slab (12 mm, 60% of thetube diameter) within the 3D volume transecting the center of the tubephantom. The data therefore represent the through slab average of thevelocities which are somewhat lower than the expected peak velocities(approx. 16 cm/s for a flow rate of 25 ml/sec).

[0057] The effects of gradient field non-uniformity on apparent flowvelocity profiles can clearly be identified and are most significant forvelocity components v_(z)(x,z) along the principal flow direction(z-direction, FIG. 2A). Note the falloff of uncorrected velocities atthe edges of the FOV (arrows), which are reduced with the generalized PCreconstruction, as demonstrated by the high intensities near the edgesin the difference image. This behavior is also consistent with thespatial variation of relative gradient strength of the respectivegradient axis (z) used for velocity encoding as depicted in FIG. 1A fora similar coronal location (top right, d_(z)(x,z) for y=0 mm).

[0058] However, apparent velocities along directions orthogonal to thetube axis such as v_(y)(x,z) (y-direction, FIG. 2B) are also affected bygradient field non-uniformity. The non-zero apparent velocities along yfor the standard reconstruction (FIG. 2B, left plot) can be explained bythe direction error associated with the gradient field used for encodingof the y-velocities. Any deviation from the nominal direction willresult in an unwanted contribution of the strong z-velocities to theintended y-encoding direction. The increase of v_(y) with distance fromthe iso-center agrees well with the spatial dependence of the of angulardeviation from the nominal encoding direction depicted in FIG. 1C for asimilar coronal location (top, mid, θ_(y)(x,z) for y=0 mm).

[0059] Both figures demonstrate that the gradient field induced errorscan be corrected with the proposed generalized reconstruction (equation7). Relatively constant laminar velocity profiles along z, as expectedfrom steady flow conditions and phantom location, can be regained alongthe entire tube after employing the generalized phase contrastreconstruction. The remaining velocity gradient for v_(y) may be due toresidual deviations between the model and real gradient fields.

[0060] In addition, it is evident that the divergence from nominalvelocity increases with distance from the iso-center of the magnet anddeviations can be as large as 25% for this phantom location.

[0061] The effects of standard and generalized reconstruction arefurther illustrated in the graphs in FIGS. 3A-3C for v_(x), v_(y), andv_(z). Velocities for all encoding directions (x,y and z) are shownthrough the center of the tube as indicated by the dashed line in FIG.2B (right). In FIGS. 3A-3C, 3D phase contrast MRI with three-directionalvelocity encoding includes measured velocities v_(x), v_(y) and v_(z)along cross sections through the center of the tube phantom (asindicated by the dashed line in the right plot in FIG. 2B). Each graphshows velocities for standard (solid gray line, standard PC) andgeneralized phase contrast reconstruction (thick black dashed line,generalized PC) as well as the difference between both methods (blackdot-dashed line).

Approximate Reconstruction

[0062] To estimate the error associated with the approximatereconstruction (equation (12)) compared to the generalized method(equation (7)), relative deviations from true velocities forone-directional velocity encoding were generated from simulations basedon the gradient field model. Table 2 lists the resulting maximumreconstruction errors (percentage deviation of the encoded velocitycomponent from the actual velocity component) for encoding along the zand x within two different volumes ((200 mm)³ and (400 mm)³) centered atthe iso-center of the magnet. Errors were calculated for identicalvelocity components (v_(x)=v_(y)=v_(z)) and a more practical examplecorresponding to flow velocities directed mostly along the principalencoding axis (v_(x)=v_(y)=¼v_(z) for encoding along z,v_(z)=v_(y)=¼v_(x) for encoding along x). For comparison, the errors forthe standard phase contrast processing (equation (8)) are also listed.TABLE 2 Comparison of maximum errors in encoded velocities for standardand approximate reconstruction within measurement volumes of (200 mm)³and (400 mm)³, for two cases: (1) equal velocities in all directions and(2) velocity in the encoding direction is 4 times larger than in theother two directions. Max. reconstruction error, Max. reconstructionerror, encoding using z-gradient coil^(#) encoding using x-gradientcoil^(##) (1) v_(x) = v_(y) = v_(z) (2) v_(x) = v_(y) = ¼v_(z) (1) v_(x)= v_(y) = v_(z) (2) v_(z) = v_(y) = ¼v_(x) standard approx. standardapprox. standard approx. standard approx. PC Recon PC Recon PC Recon PCRecon (200 mm)³ 5% 3% 4% 1% 31% 27% 13% 7% volume (400 mm)³ 94%  41% 72%  10%  153%  108%  70% 27%  volume

[0063] From Table 2, it is evident that standard phase contrastreconstruction for one-directional velocity encoding results insignificant errors in the encoded velocities which can be as large as150%. Especially for the z-gradient coil those errors can beconsiderably reduced using the approximate reconstruction. Even for alarge volume of (400 mm)³ most of the deviation of gradient strength isalready determined by the self terms and can therefore be correctedusing equation (12). Velocities measured using the x-gradient coil, onthe other hand, do not benefit as much from the approximate correction.Here, cross terms play a more significant role which is alsodemonstrated in the field plots for x- and y-gradient coils in FIG. 1 (xand y coils are identical but rotated). Nevertheless, if the flowdirection is located mostly along the encoding direction, even for the xor y gradient coil the maximum encoding error can be reduced byapproximately 50% and it is worthwhile performing the correction if onlyone-directional velocity data is available.

[0064] In general, the improvement in accuracy can be maximized if theencoding direction is adjusted to match the principal flow direction(e.g., for v_(x)>>v_(z), v_(y) for encoding along x).

[0065] Results of the 2D phantom experiments with one-directionalvelocity encoding are shown in FIG. 4. The theoretically predictedinhomogeneity of z-gradient fields used for the approximate description(equation (12)) of errors in the encoded velocity are shown in the leftcolumn. The surface plots of Λ_(zz)(x,z) demonstrate the relativedeviations from nominal z-gradient strength in the coronal plane forspecific tube y-positions. FIGS. 4A and 4B illustrate 2D phase contrastMRI with one-directional velocity encoding: demonstration of the effectof gradient field inhomogeneities on velocity encoding for locations ofthe tube phantom as defined by the (x,y) coordinates. Simulated andnormalized z-gradient fields (left) depict the spatial variation of therelative deviation λ_(zz)(x, z) from the ideal gradient strength incoronal planes at the respective phantom locations for a FOV of 400 mm.Measured through plane velocities are compared to predicted measurementerrors (dark lines) for both phantom locations. Note that the slicelocations for the measured velocities have been corrected for gradientfield distortions.

[0066] Mean measured velocities {overscore (v)}_(z,standard PC) werescaled to represent local deviations from ideal gradients (i.e., forΛ_(ii)=1) and could thus be directly compared to relative changes ingradient strength. Results are displayed as a function of slice location(z) (right column). The thick solid lines in all graphs indicate thez-gradient according to the model for corresponding phantom locations(given as (x,y) in mm).

[0067] Since steady flow was used for all experiments, the true meanthrough plane velocities are expected to be constant as a function ofspatial location (z) along the tube. As shown in FIG. 4, the deviationsof the velocities measured demonstrates the effect of gradient fieldnon-uniformity and corresponded well to the theoretically predictederrors for different phantom locations. This is further demonstrated inFIG. 5 which shows excellent correlation between predicted velocityerrors and the measured through plane velocities (linear regression,slope=0.92+/−0.06, intercept=0.05+/−0.51, r=0.94). In FIG. 5, 2D phasecontrast MRI with one-directional velocity encoding includes correlationof measured versus model based deviations Δ{overscore (v)}_(z) fromthrough plane velocities. Velocities were scaled to represent percentagedeviations from conditions with ideal gradients. The solid linerepresents the result of a linear regression.

[0068] Since the measured errors are in agreement with predictions, anapproximate error correction can be performed by a simple scalarmultiplication of the measured velocities by spatially dependent valuesΛ_(zz)(z) derived from the theoretical model of the gradient fields.

[0069] With the demand for high gradient strength and slew rates,systems with short gradient coils are being used and gradientnon-uniformity can be a significant problem. Corrections are availableto remove spatial distortions that result from imperfect spatialencoding. However, errors introduced into phase contrast measurementshave not been thoroughly reported.

[0070] The effect of gradient field imperfections on strength andorientation of velocity encoding demonstrates the necessity for a moregeneralized reconstruction of PC-MRI data. In general, accuratereconstruction of magnitude and direction of velocities can be achievedif three-directional velocity encoding is performed and an appropriatecorrection that includes the direction error is applied.

[0071] The simulations and velocity measurements demonstrate thatgradient field inhomogeneities can introduce significant errors inmagnitude and direction in quantitative phase contrast flow or velocitymeasurements. Comparison of standard and generalized reconstruction ofthree-directional phase contrast data demonstrate that a gradient fieldmodel based on a spherical harmonic expansion can be used to compensatefor magnitude and angular deviations in first gradient moments and thusvelocity encoding. The relevant parameters for the field model aresystem dependent.

[0072] From both model and measurements it is evident that the errors invelocity encoding strength increase with distance from the iso-center ofthe magnet while the angular deviation demonstrates more complexpatterns. In general, through plane velocity measurements with a singleslice placed at or near the center of the magnet and with the vessel ofinterest close to the center of the FOV are largely insensitive toencoding errors. The situation differs considerably, however, ifmultiple slices are acquired within a single acquisition or if bloodflow is to be analyzed within a larger plane or volume.

[0073] The effect on the magnitude of encoded velocities stronglydepends on the off center location of the object and can, in our system,be as large as 10% within a (200 mm)³ volume and 65% within a (400 mm)³volume. Angular deviations, can result in encoding direction errors upto 16° within a (200 mm)³ volume and 45° within a (400 mm)³ volume.

[0074] Especially for measurements within smaller volumes, the directionerror (contribution from cross terms) may be ignored and encodedvelocities can be corrected by a simplified method (equation (12)). Thissimplified version effectively ignores any deviation from the originalencoding direction but has the advantage that one-directional velocitymeasurements can be corrected since cross terms are neglected.

[0075] Although accurate reconstruction would require the fullthree-directional velocity information the second set of phantomexperiments demonstrate that the error in through plane velocities couldbe successfully described by this approximate correction. Even thoughnot all gradient field deviations are accounted for and some encodingerror remains, significant improvement could be achieved if compared tothe uncorrected data.

[0076] It has to be noted, however, that the excellent correlationobtained in the analysis of the phantom data is also due to the specialgeometry used herein. Since, according to the gradient field model, theencoded velocity is not in general parallel to the true velocity,directional errors and thus cross-terms may be more significant fordifferent geometries with major flow components not along the principalvelocity encoding direction.

[0077] Gradient non-uniformities also cause distortion of the excitedslice in 2D imaging (gradient non-linearities causes an RF pulse toexcite spins that are not within a flat plane, ‘potato-chip shape’).

[0078] The local normal to the excited slice is parallel to the localgradient direction. Thus, if flow encoding is in the “slice direction”,the flow encoding is parallel to the local normal to the slice, as wouldbe desired. The slice distortion could cause the image to not containthe desired anatomy, but the effect is difficult to estimate since theseerrors are dependent on the local gradient field deviations and as suchhighly dependent on the location of the anatomy within the slice, aswell as the slice plane location itself. However, assuming the imagecontains the vessel of interest, if the measured slice is corrected forspatial distortion and the velocity measurement is corrected for theerror in the flow encoding strength, the measured volume flow rateshould be correct. Note that in this case, correction for the error inthe flow encoding direction would be counterproductive.

[0079] In summary, deviations related to gradient field inhomogeneitiesneed to be taken into account if such gradient fields are used to encodeflow as in phase contrast MRI. Significant improvement in velocityquantification can be achieved by using the known field distortions tocorrect the measured phase shifts. After spatial unwarping, errors inquantitative phase contrast measurements can be greatly reduced withrelatively little effort by using the generalized reconstructionalgorithm in accordance with the invention.

[0080] While the invention has been described with reference to specificembodiments, the description is illustrative of the invention and is notto be construed as limiting the invention. Various modifications andapplications may occur to those skilled in the art without departingfrom the true scope and spirit of the invention as defined by theappended claims.

Appendix

[0081] Gradient field model:

[0082] To model the spatial dependence of the field generated by eachgradient coil (i=x, y, or z coil) a spherical harmonic expansion isemployed: $\begin{matrix}{{B_{z,{model}}^{(i)} = {(r) = {{\sum\limits_{k = 0}^{\infty}{\sum\limits_{l = 0}^{k}{\alpha_{kl}^{i}{H_{kl}^{C}(r)}}}} + {\beta_{kl}^{i}{H_{kl}^{S}(r)}}}}},{i = x},{y\quad {and}\quad {z.}}} & ({A1})\end{matrix}$

[0083] H_(kl) ^(C)(r) and H_(kl) ^(S)(r) are the solid sphericalharmonic basis functions and α,β represent coefficients that arespecific for an individual MR unit and typically provided by themanufacturer. For each gradient coil three terms are included in theexpansion. From the model, it can be shown (see FIG. 1) that anenergized gradient coil produces not only changes in B_(z) along thenominal gradient orientation but also in the two orthogonal directions.

Generalize Reconstruction—Common Special Cases

[0084] In this section examples of the application of the generalizedreconstruction (equation (7)) for typical velocity encoding in phasecontrast MRI are given.

[0085] If one-directional velocity encoding is performed with anencoding gradient G_(z,ideal)(τ) along the nominal z-axis only(G_(x,ideal)(τ)=G_(y,ideal)(τ)=0) the true first gradient moment isgiven by: $\begin{matrix}\begin{matrix}{{M_{1}(r)} = \begin{pmatrix}{M_{1,x}(r)} \\{M_{1,y}(r)} \\{M_{1,z}(r)}\end{pmatrix}} \\{= {\int_{0}^{TE}{\begin{pmatrix}{\Lambda_{xx}(r)} & {\Lambda_{xy}(r)} & {\Lambda_{xz}(r)} \\{\Lambda_{yx}(r)} & {\Lambda_{yy}(r)} & {\Lambda_{yz}(r)} \\{\Lambda_{zx}(r)} & {\Lambda_{zy}(r)} & {\Lambda_{zy}(r)}\end{pmatrix}\begin{pmatrix}0 \\0 \\{G_{z,{ideal}}(\tau)}\end{pmatrix}\tau {\tau}}}} \\{= {{\begin{pmatrix}{\Lambda_{xz}(r)} \\{\Lambda_{yz}(r)} \\{\Lambda_{zz}(r)}\end{pmatrix}{\int_{0}^{TE}{{G_{z,{ideal}}(\tau)}\tau {\tau}}}} = {\begin{pmatrix}{\Lambda_{xz}(r)} \\{\Lambda_{yz}(r)} \\{\Lambda_{zz}(r)}\end{pmatrix}{M_{1,z,{ideal}}(\tau)}}}}\end{matrix} & ({A2})\end{matrix}$

[0086] As a result, the actually measured phase also includes cross-terncontributions from orthogonal gradient terms determined by Λ_(zx)(r) andΛ_(zy)(r). The measured phase difference Δφ_(z)(r) in this case is thusgiven by $\begin{matrix}\begin{matrix}{{{\Delta\varphi}_{z}(r)} = {{v(r)}^{T}{\gamma \begin{pmatrix}{\Lambda_{xz}(r)} \\{\Lambda_{yz}(r)} \\{\Lambda_{zz}(r)}\end{pmatrix}}\Delta \quad M_{1,z,{ideal}}}} \\{= {{\gamma \left\lbrack {{{v_{x}(r)}{\Lambda_{xz}(r)}} + {{v_{y}(r)}{\Lambda_{yz}(r)}} + {{v_{z}(r)}{\Lambda_{zz}(r)}}} \right\rbrack}\Delta \quad M_{1,z,{ideal}}}}\end{matrix} & ({A3})\end{matrix}$

[0087] Therefore, without additional information (e.g., it is known thatv_(x)(r) and v_(y)(r) are negligible) the actual velocity cannot berecovered from the data.

[0088] For a typical phase contrast experiment with velocity encodingalong the three orthogonal gradient coil axes x, y and z (N=3) theencoding matrix reduces to a (3×3) diagonal matrix. $\begin{matrix}{\Omega = \begin{pmatrix}{\Delta \quad M_{1,x,{ideal}}^{(1)}} & 0 & 0 \\0 & {\Delta \quad M_{1,y,{ideal}}^{(2)}} & 0 \\0 & 0 & {\Delta \quad M_{1,z,{ideal}}^{(3)}}\end{pmatrix}} & ({A4})\end{matrix}$

[0089] For identical velocity sensitivities ΔM=ΔM_(1,x,ideal)⁽¹⁾=ΔM_(1,y,ideal) ⁽²⁾=ΔM_(1,z,ideal) ⁽³⁾ along all encoding axes themeasured phase differences are given by ΔΦ(r)=γΔM[Λ(r)]^(T)v(r) and thecalculation of the true velocities (equation (7)) can be simplified to$\begin{matrix}{{v(r)} = {\frac{1}{{\gamma\Delta}\quad M}\begin{matrix}\left\lbrack {\Lambda (r)} \right\rbrack^{- T}\end{matrix}{{{\Delta\Phi}(r)}.}}} & ({A5})\end{matrix}$

1. A method for measuring velocities in phase contrast magneticresonance imaging in the presence of magnetic field gradient fieldinhomogeneities comprising the steps of: a) encoding velocityinformation in at least three different encoding directions andmeasuring phase differences (ΔΦ(r)) in each direction, b) derivingactual first moments along the encoding directions using a magneticfield deviation matrix (Λ(r)), c) providing an (N×3) encoding matrix, Ω,where N is the number of encoding directions, based on the ideal firstmoments along the N directions, and d) computing true velocities as afunction of the gyromagnetic ratio, γ, magnetic field deviation matrix(Λ(r)), encoding matrix, Ω, and a phase difference vector (ΔΦ(r)). 2.The method as defined by claim 1 wherein true velocity, v(r), iscomputed as follows: ${v(r)} = {\frac{1}{\gamma}\begin{matrix}\left\lbrack {{\Lambda (r)}\Omega} \right\rbrack^{- T}\end{matrix}{{{\Delta\Phi}(r)}.}}$


3. The method as defined by claim 2 wherein the magnetic field gradientsare produced using gradient coils, and the magnetic field deviationmatrix is given by ${\Lambda (r)} = \begin{pmatrix}{\Lambda_{xx}(r)} & {\Lambda_{xy}(r)} & {\Lambda_{xz}(r)} \\{\Lambda_{yx}(r)} & {\Lambda_{yy}(r)} & {\Lambda_{yz}(r)} \\{\Lambda_{zx}(r)} & {\Lambda_{zy}(r)} & {\Lambda_{zy}(r)}\end{pmatrix}$

where the lambda terms are the deviations of the magnetic field producedby each gradient coil from the corresponding ideal gradient.
 4. Themethod as defined by claim 2 wherein the encoding matrix, Ω, is givenby: Ω=(ΔM _(1,ideal) ⁽¹⁾ , ΔM _(1,ideal) ⁽²⁾ , . . . , ΔM _(1,ideal)^((N))).
 5. The method as defined by claim 4 wherein velocitysensitivities ΔM are identical along all encoding directions and thetrue velocity is given by:${v(r)} = {\frac{1}{{\gamma\Delta}\quad M}\begin{matrix}\left\lbrack {\Lambda (r)} \right\rbrack^{- T}\end{matrix}{{{\Delta\Phi}(r)}.}}$


6. A method for computing velocities (v_(⊥)(r)) in phase contrastmagnetic resonance imaging in the presence of magnetic gradient fieldinhomogeneities comprising the steps of: a) measuring velocityinformation including phase difference Δφ₁(r) along an encodingdirection and the absolute gradient deviation ∥Λ(r)n∥ in that direction,b) calculating flow encoding strength,$\Delta \quad M_{1,{ideal}}^{(1)}\sqrt{\sum\limits_{j = 1}^{3}\left\lbrack {\sum\limits_{i = 1}^{3}{n_{i}{\Lambda_{ji}(r)}}} \right\rbrack^{2}}$

based on relative magnetic field gradient deviations and first momentdifference, ΔM_(1,ideal) ⁽¹⁾, and c) computing the velocity usinggyromagnetic ratio, γ, as:${v_{\bot}(r)} = {\frac{{\Delta\varphi}_{1}(r)}{{\gamma\Delta}\quad M_{1,{ideal}}\sqrt{\sum\limits_{j = 1}^{3}\left\lbrack {\sum\limits_{i = 1}^{3}{n_{i}{\Lambda_{ji}(r)}}} \right\rbrack^{2}}}.}$